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SUMMARY 

Two computer programs coded in FORTRAN IV are described for plane collision- 
less positive-ion and electron emission sheaths. Given the emitter temperature, emit- 
ter work function, atomic ionization potential, plasma electron and ion number density, 
and plasma electron, ion, and atom temperatures, the programs compute current den- 
sities, potential drop through the sheath, charge density, electron field, and sheath dis- 
tance. 


INTRODUCTION 

Methods to calculate the properties of sheaths between a thermally ionized plasma 
and a field- modified emitter of electrons, ions, and atoms are programmed in 
FORTRAN IV for the IBM 7094. References 1 and 2 give the background, theory, and 
many results for emission sheaths and also graphic correlations for sheath character- 
istics in cesium plasmas. The examples herein are specific for a cesium plasma. The 
programs are general and convert readily to other plasmas by replacing the minimum 
mean free path, vapor pressure, ionization potential, and masses of the atom and ion of 
cesium with those of the desired chemical. Computing procedures for the positive-ion 
and electron sheaths appear as FORTRAN programs in appendixes A and B. Appendix C 
defines the FORTRAN variables and corresponding output labels. Appendix D defines 
additional FORTRAN variables and constants. Flow diagrams of the programs appear in 
figures 1 and 2. 


SYMBOLS 


E electric field 

e electronic unit charge 

h Planck' s constant 

I ionization potential 

J overall net current 

j current density or net particle current density 

m particle mass 

N particle number density 

p pressure 

T temperature 

V potential 

AV potential relative to plasma potential 

X distance from emitter 

k Boltzmann constant 

X mean free path for charge exchange in emitter temperature range 

X D plasma Debye length, =6. 9 ( T e p/ N ep ) 1//2 
X DE emission Debye length, ~6. 0 (T E /N e p) 1//2 
<p work function 

<P Q plasma potential for near-equilibrium positive-ion sheath 

Subscripts: 
a atom 

E emitter 

e electron 

i ion 

p plasma 

S overall sheath 

AV potential relative to plasma potential 

vp vapor pressure 
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DESCRIPTION OF INPUT 


The input for one case requires one value each of T F> <P, N pn> T mV and T in‘ The 

N ^ q 

ranges of the conditions used in the calculations for cesium plasmas are 1400 to 
2400° K for the electrode temperature, 1600° to 2400° K for the plasma atom and ion 
temperatures, 1700° to 2700° K for the plasma electron temperature. 10 12 to lO 1 ^ for 
electron number density, and 1. 5 to 5. 0 volts for the work function. The work function 
is an assigned variable; the equilibrium value for no sheath was obtained by solving for 
the Saha- Langmuir null point 


cp = — ( — In T - In N 


ep 


o 27rm k 
+ £ In S_ 


+ In 2 


(1) 


and choosing three or four values of cp above this point for the positive ion sheath cal- 
culations, and three or four values below the null point for the electron sheath calcula- 
tions. 

More than one case can be run during one machine access. On the 7094 II-DCS, ap- 
proximately six cases can be run in 1 minute. 

Input card preparations for the Positive-Ion and Electron Sheath Programs are 
identical. A single digit integer in column 5 of the first card is IWRITE, the variable 
that controls the desired output. If IWRITE = 0, results for Richardson- Dushman (refs. 

3 to 5) and Schottky (ref. 6) are printed. If IWRITE = 1, the Richardson- Dushman output 
is eliminated. The remaining cards are in five groups. In each group the first card (an 
integer in cc4-5) indicates the number of values to be read from the successive cards of 
that group (Format 8E10. 2): 


Group 

Fortran variable 
name of count 
(-10) 

Fortran 

array 

variable 

Description 

Equation 

variable 

I 

n 

TE 

Emitter 

temperature 

t e 

II 

jj 

PHI 

Emitter work 
function 

ecp 

in 

KK 

EPN 

Plasma electron 
number density 

N ep 

IV 

LL 

TEP 

Plasma electron 
temperature 

T ep 

V 

LA 

TIPP 

Plasma ion 
temperature 

T 1P 
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A second set of five groups of input data cards may follow; IWRITE is not included after 
the first set is read in. 

In order to convert the programs from a cesium plasma to a plasma of another ele- 
ment, replace the ionization potential, AI, (FORTRAN statement 301; mass of the atom, 
AM, (302); vapor pressure, PTEST, (303); and minimum mean free path, AMTEST, 
(304) with those of the desired element. 


METHOD OF CALCULATION 


The programs described deal specifically with a cesium plasma; therefore, the 
tests and constants used are characteristic of cesium. As mentioned earlier in the sec- 
tion DESCRIPTION OF INPUT, the programs can be converted to handle another chemical 

if desired. The calculations in the programs assume m. = m„ and N. = N__. 

i ct ip ep 

Before major calculations begin, the input values must satisfy two tests. First, the 
Debye shielding length (AMDA) must be less than the low-energy mean free path for ce- 
sium charge exchange (AMTEST). This enables the model to be collisionless. Second, 
the vapor pressure of cesium at the emitter temperature (PTEST) must be less than the 
effective plasma pressure (PPT). The equations for these variables are as follows: 


AMDA = 


(4rr)(0. 511xl0 6 )(2. 82 x10" 13 ) 


1/2 At \l/2 
ep 


cm 


i N 


ep> 


( 2 ) 


AMTEST = 10 12 N^J, cm 


(3) 


PTEST = 


anti log ( - 392 °- 38 - 0. 51978 log Tg + 10. 71914' 
T E 


133. 322 


-, torr 


(4) 


pp T _ (82. 06)(760) / 2N T N T \ torr 
9 o v ep ip ap ap ; ’ 

6. 023x10^ 


(5) 


(133. 322 (N/m 2 )/torr) 

If the input for one case passes the preceding tests, calculations continue; otherwise, 
another case with a new N g p is begun. 


4 


The Richardson-Dushman equation and e<p' = e<p yield the current densities if the 
effects of the sheath field at the emitter, E £ are omitted. 

Electron emission current density: 


J eE 


120 T 


E 




( 6 ) 


Plasma current density: 
Electron: 


j 


ep 


1. 602x10' 19 N ep 
2 



Ion: 


J ip 


1. 602x10' 19 N. 
2 



Atom: 


1. 602x10' 19 Np I 


'2kT 


. 1/2 


ap 




J ap 


27rm e KT ep 


1 3/2 


exp 1 - 


el 


K T 


epy 


(7) 


(8) 


(9) 


Equations (6) to (9) are the same for positive-ion and electron sheaths. 

In order to calculate the current densities through the sheath, a need arises for the 
overall potential sheath drop AV g . A first approximation is made by using the charge 
density equations P^y and conditions at the plasma, sheath interface, where AV = 0. 
The following equations define the charge density, equation (10) for the positive- ion 
sheath and equation (11) for the electron sheath. 


5 
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Electron sheath (AVg negative): 


J eE 


AV 


1 + erf ( - 


e A V 
kT t 


a/2 


exp 


e(AVg - AV)' 


kT 


'2kT t 


a/2 


E 


7rm, 


J ep 


1 - erf I 


e AV 


a/2 


kT 


e P/ 


exp I- 


e AV 


kT 


epy 


' 2KT ep' 1/2 
7rrn 



r 

e(AV q - AV)' 

1 / 2 ) 



1 - erf 

kT e 

I 

► exp - 


v. 


J 

*- 


e(AV g - AV) 


kT, 


(2kT t 


a/2 


jrnr 


J ip 


r 

1 + erf 

e(AV g - AV)' 

1/2| 

/cT.,, 

( 

l 

ip 

J 


►exp 


e AV 


k T. 




(2kT- 
* ± ip 


a/2 




When AV = 0, equations (10) and (11) reduce to the charge density in the plasma p 



Positive-ion sheath: 



Electron sheath; 



V m i/ 


In equation (12), the third term equals one-half the ion charge density N. ; equation (16) 

■*-ir 

is substituted for In equation (13), the first term equals one-half the electron 

charge density 1ST; therefore, an approximation of AV C can be made. 

^ 6p o 
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Positive -ion sheath: 



Electron sheath; 


AV S 


kT, 


1. 602x10" 19 N 


ep 


In 


2j 


eE 



(15) 


A test is made on AVg to assure that it is positive for the ion sheath and negative for 
the electron sheath. If the test is not satisfied, a new case is begun by using another 
value of cp . 

With the first estimate of AVg, one can solve for the current densities from the 
emitter and the net negative current densities through the sheath into the plasma. 

Ion current density from emitter: 

Positive-ion sheath: 


J iE “ 


J ip + ] ap 


2 exp 


e (I - <?' ) 


K T 


+ exp I - 


E 


e AV C 
kT t 


(16) 


Electron sheath: 


•*ap + iip ex P I 


^iE 


fe AVr 


KT ip 


1 + 2 exp 


e(I - <P') 

K 


(17) 


9 



Atom current density from emitter; 
Positive- ion sheath: 


J aE 


•’ip + •’ap 


1 + - exp 
2 


e(AV g - <p' + I) 


k T t 


Electron sheath: 


e AV C 


V + lip exp i 


K T- 


J aE 


iP 


1 + — exp 
2 


_ e(I - <p») 

kT e 


Net negative current density through sheath into plasma: 
Electron: 

Positive-ion sheath: 


e AV C 


)e=leE~ Jep ex P|- 


kT 


ep y 


Electron sheath: 


3 e = i e E exp 


e AVg t 

kT t l' 3ep 


Ion: 


Positive-ion sheath: 


h = hp~hE exp l“ 


e AV £ 

kTtp 


( 18 ) 


( 19 ) 


( 20 ) 


( 21 ) 


( 22 ) 
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Electron sheath: 


h = itp ex P 


'e AV, 


- 3 


iE 


(23) 


Atom: 

Positive-ion and electron sheath: 


^a ^aE " ^ap (^) 

Overall: 

J = 3 e + jj (25) 

The overall potential drop AVg is then divided into 20 equal increments AV, and 
the net negative charge density in the sheath equations (10) or (11) is calculated for in- 
creasing potentials. The error function in equations (10) and (11) was evaluated by using 
the Lewis library subroutine for error functions. 

Integration of the net negative charge density yields the electron field in the sheath 
E^y. The calculations use the sheath fields and voltages given in terms of electron po- 
tential; therefore, positive signs correspond to ion sheath and negative signs correspond 
to electron sheath: 



(28) 


The numerical integration of p^y is performed by using the trapezoidal rule based on 
20 equal increments of AV from 0. 0 to AVg. 

With the field at the Schottky emitter E^ = E^y , a Schottky correction 

e(p' = e<p - e|(0. 511xl0^)(2. 82x10”^ Eg)J^^ usec * * n e Q. ua tion (6) for the ion sheath 

program, and ecp' = ecp + e|(-0. 511xl0^)(2. 82x10”^ Eg)J^ is substituted in equa- 
tion (17) for the electron sheath program. The Schottky correction is used in the recal- 
culation of the values of current densities, overall potential drop, charge density, and 
electron field. This cycle continues until no boundary current density affected by the 
Schottky correction (for example CAT in the Positive-Ion Sheath Program) changes by 
more than 0. 1 percent of the smallest boundary (plasma or emitter) current density 
(TENT). Once the preceding test is satisfied, the sheath distances X are computed: 
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dAV 

E 


(27) 


^AV 


/O /*0 

: c */av 


As in the electron field calculations, the trapezoidal rule is used for numerical integra- 
tion. However, at zero AV, E equals zero, and equation (27) is not defined. The pro- 
grams deal with this singularity by using E/2 at AVg/20 as the average for the first 
increment and compute a finite X at this point. 


DESCRIPTION OF OUTPUT 

Examples of the output are shown in figure 3 for the positive-ion sheath and in fig- 
ure 4 for the electron sheath. Following the numerical values (figs. 3(a) and 4(a)) are 
plots which show sheath voltages, fields, and charge densities as functions of distances 
from the electrode. The sheath properties vary with the plasma electron concentration 
N g p, electron and ion temperatures T g p and T^, electrode work function <p, and 
emitter temperature T £ . The first set of five plots (plotted by the method of ref. 7) 
shows the Richardson-Dushman and Schottky results (figs. 3(b) and 4(b)); the second set 
shows only the Schottky results (figs. 3(c) and 4(c)). If IWRITE = 1, the second set of 
five plots is eliminated. The sample output corresponds to the input listed at the end of 
the Fortran programs (appendixes A and B). 

An explanation of the output labels and corresponding FORTRAN variables is given 
in appendix C. Appendix D gives an alphabetical listing of important FORTRAN vari- 
ables, not defined in appendix C. The recorded plasma and electrode properties, over- 
all sheath characteristics, and tabulated and graphic incremental data are sufficient to 
detail the structure and transport of an emission sheath. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, December 4, 1967, 

120-33-02-01-22. 
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APPENDIX A 


POSITIVE -ION SHEATH PROGRAM LISTING 


$ I BFTC ION 
C 

C A PLANAR ION SHEATH BETWEEN AN EMITTER AND A 

C NEAR - EQUILIBRIUM PLASMA 

C 

COMMON /MV/ PI, EM, AK, AM 

COMMON /MP/ XDVI 50) , XDS ( 50 ) , RHODOP ( 50 ) , RHDS ( 50 ) , EDV ( 50 ) , EDS ( 50 ) , 
1DV( 50) ,OVSAVE( 50) ,RHOEOP(50) , RHES ( 50 ) , RHO I OP < 50 ) , RHI S ( 50 ) , NBB 
2, I WRITE 

DIMENSION TE(10),PHI ( 10 ) , EPN { 10 ) , TEP ( 10) ,EJ(30) ,EPJ(30) ,PIJ(30), 
1APJI30) ,PP(30) ,DVS( 100) ,E1J(30) ,AEJ(30) ,CJ(30) ,RH0D(40) , 
2ESAVE(40 ) ,XSAVE (40) ,RH0E(40) , T I PP ( 1 0 ) , E JS ( 50 ) , AMDAI 10) ,SJE( 20) , 
3SJA(20),SDJA(20) , SDJE < 20 ) , SDJ I ( 20 ) , J 1 (100 ) , J2 ( 100 ) , 

4J3( 100) , J4< 100) ,RHOI (40) , AMDATEI 10) , Y < 40 ) , SJ I ( 20 ) 

REAL J1 , J2 ,J3,J4, JA, JB, JC 
C 

C A I , IONIZATION POTENTIAL 

C 

301 A I = 3.893 

AK = 8.617E-5 
AKT = 2.0*AK 

EM = .511E+6/I <2.998E+10)**2) 

PI = 3.14159 
P I EM = P I *EM 
SE = 1.602E-19 
C 

C MA , MASS OE ATOM 

C 

302 AM = 931.478E+6*132.9/((2.998E+10)**2) 

PI AM = PI 4AM 

API = 6.6256E-34/SE 
AP2 = AM 

Cl = SORTIAK/ (4.04PI4.511E+642.82E-13) ) 

C2 = AP 1442 

C3 = 82. 064760. /6.023E+23 
C4 = .511E+642.82E-13 
C5 = -72.04PI41 .OE+11 
C 

C IWRITE = 1 FOR SCHOTTKY OUTPUT 

C IWRITE = 0 FOR RICHARDSON - DUSHMAN AND SCHOTTKY OUTPUT 

C 

READ (5,12) IWRITE 
C 

C READ INPUT TE , PHI , EPN , TEP , TIP 

C 

1 READ (5,2) I I , (TE< I ) , 1=1, I I ) 

READ (5,2) JJ, (PHI (J) , J = 1 , J J ) 

READ (5,2) KK, (EPN(K) ,K=1,KK) 

READ (5,2) LL, (TEP(L) »L=1*LL) 

READ (5,2) LA, (TIPP(LB) ,LB=1,LA) 

DO 1000 I =1,11 
C 

C COMPUTE PTEST , VAPOR PRESSURE 

C 

303 PTEST = 10.**(-3920.38/TE( I ) -. 5 1 978 1 4AL0G1 0 ( TE ( I ) )+10. 71914)/ 

1 133.322 

AKTE = AK4TE ( I ) 
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TAKTE = 2 . O^AKTE 

SA = SORT ( T AKT t /P I AM ) 

FD = SA**2 

CI1 = P I EM*TAKT b 

C I 2 = SORT (TAKTE/PIEM) 

TES = T E ( I ) **2 
DO 995 LB = 1 * L A 
TIP = TIPP(LB) 

CL B 1 = SORT (AKT*TIP/PIAM) 

SET TAP = TIP FCJR THESE SOLUTIONS 

TAP = TIP 
DO 990 L = 1 , L L 
AKTP = AK*TEP(L) 

TAKTP = 2.0-AKTP 
SB = TAKTP*PIEM 
SB 1 5 = SB**1.5 
CL 1 = EXP ( A I /AKTP ) 

CL2 = SORT (TAKTP/PIEM) 

FB = CL2**2 
DO 980 K=1,KK 

COMPUTE PHIZZ , AMDA , AMDATE , APN , AMTEST , PPT 
APN = (EPN(K>*C2>*(EPN(K)*(AP1) )*CL1/SB15 

CK1 = AKTE*ALOG< .5*( ( T I P/TE ( I ) ) **. 5+ ( APN/EPNI K ) )* ( TAP/ TE < I ) )**. 5 
1 - 1 . 0 ) ) 

SEE = SE*EPN(K) 

PHIZZ = -AK.TE*ALOG ( (SEE/(240.*TES) ) *C I 2 ) 

AMDA(K) = C1*SQRT(TEP(L)/EPN(K) ) 

AMDATE(K) = C 1 *SORT ( T E ( I) / E PN ( K ) ) 

AMTEST, MINIMUM MEAN FREE PATH 

304 AMTEST = 1.0E+12/APN 

PPT = C3*(2.0*EPN(K)*TIP+APN*TAP) 

TEST ON AMDA , AND AMTEST 

IF ( AMDA (K).LT. AMTEST) GO TO 9 

WRITE (6,7) AMDA(K), AMTEST, TEII ),EPN(K), TEPID 

GO TO 980 

TEST ON PPT AND PTEST 
9 IF (PPT -PTEST) 40,40,5 

5 WRITE (6,6) PPT , PTEST , TE ( I ) , EPN ( K ) , TEP ( L ) 

GO TO 980 
40 CONTINUE 

DO 970 J = 1,JJ 
KOUNT = 0 
METS = 0 
IPRINT = IWRITE 
3 PHD = PHI ( J) 

COMPUTE EJ , EPJ , PIJ , APJ , PP 

EJ ( J ) = 120.*TES*EXP (-PHD/AKTE ) 

EPJ(J) = SEE*.5*CL2 
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PIJ(J) = « 5*S EE^CLB 1 

APJ(J) = (SE*EPN(K)**2*AP1 )*(CL1*C2 )*SQRT ( AKT #fAP/(PI*AP2)) 
1/<S615*2.0) 

PP ( J ) = C3*{ EPN(K)*(TEP( L)+TIP*1.0+TAP*( ( EPN ( K ) *C2 ) *CL 1 *A P 1 ) / SB 1 5 ) 
1) 

IF { IWRITE.EQ.l.OR.METS.EQ.l) GU TU 38 
WRITE (6,799) 

8 WRITE (6,800 ) A I , TE ( I ) , PH I ( J ) , EPN ( K ) , TEP ( L ) , T I P , AMDA ( K ) , PTE ST 
1 , AMDATE ( K ) 

38 CONTINUE 

KOUNT = KOUNT + 1 

COMPUTE DVS 

106 DVS(J) = CK1+PHD-AI 

IF (METS.EO.O) DVSRD = OVS(J) 

PHIZ = PHI ( J) - DVSRD 
OROK = ABS ( DVSRD/ AKT P ) 

DVSPZ = DVS ( J ) / DVSRD 
115 IF ( DVS ( J )) 108,106,102 
108 WRITE (6,109) 

GO TO 970 
106 WRITE (6,107) 

GO TO 970 

COMPUTE E I J , AEJ , SJE , SJI , CJ , SDJA , SDJE , SDJI , J1 
J2 , J3 , J4, SJA 

102 EIJ(J) = <PIJ< J)+APJ( J) )/<2.0*EXP( (AI-PHD)/AKTE)+EXP(-DVS( J)/AKTE) 
1 ) 

A E J ( J ) = ( P I J ( J) + APJ(J ) )/( 1.0+.5*EXP(-( DVS< J)-PHD+AI ) / AKTE ) ) 

SJE(J) = EJ { J ) -EPJ ( J ) *EXP ( —DVS ( J ) / AKTP ) 

SJI(J) = -EIJ (J ) *EXP { -DVS ( J )/AKTE)+PIJ(J ) 

C J ( J ) = S J E ( J ) + S J I ( J) 

SJA(J) = AEJ( J)-APJ( J) 

SDJA(J) = $J A ( J ) /APJ ( J ) 

SDJE(J) = S J E ( J ) / EP J ( J ) 

SDJI (J) = SJI ( J ) / P I J ( J ) 

J1 (KOUNT) = E J ( J ) 

J2 ( KOUNT ) = EPJ( J)*EXP(-DVS(J)/AK1P) 

J 3 ( KOUNT ) = EIJ{ J)*EXP(-DVS(J)/AKTE) 

J 4 { KOUiNT ) = PIJ(J) 

IN = 20 

COMPUTE - IN - VALUES OF DV FROM 0.0 TO DVS 
IN = 20 

112 DVI = DVS ( J ) / FLOAT (IN) 

D V ( 1 ) = 0.0 

DO 120 NB = 1 , I N 

NBB = NB + 1 

DV(NBB) = DV(NB) + DVI 

120 CONTINUE 

COMPUTE RHOEOP , RHOIOP , RHODOP 
ERF IS THE ERROR FUNCTION SUBROUTINE 
VIN IS A SUBROUTINE TO CALCULATE SORT (2KT/PI*M)*EXP(-DV/KT)/ 

( 1 .O-ERF ( SORT ( DV/KT ) ) ) 

121 DO 200 NB = 1 , N B 8 
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X = SORT ( DV < NB >/AKTE ) 

SRC = ERF ( X ) 

X = SORT ( (DVS( J)-DV(NB) J/AKTP) 

ERO = ERF ( X ) 

DV P = DVS( J)-DV(NB) 

IF (DVP.EO.O.O) GO TO 160 
CALL VIN ( 1 , DVP , TE ( I ) , ANS) 

AVI = ANS 
GO TO 161 

160 AVI = C 1 2 

161 CONTINUE 
DVP = DV(NB) 

IF (DVP.EO.O.O) GO TO 165 
CALL VIN (2, DVP, TIP, ANS) 

AV2 = ANS 
GO TO 166 

165 AV2 = CLB1 

166 CONTINUE 

FA = EXP(-DV(NB)/AKTP) 

AE1 = EPJ( J )*FA*( 1 .O+ERD)/ ( FB**.5) 

FC = EXP(-(DVS( J l-DV(NB) l/AKTE) 

AE2 = E I J( J )*FC*( 1.0+ERC ) / ( FD*#.5> 

AE3 = E J ( J ) / A V 1 
AE4 = P I J ( J ) / AV2 
185 RHOE(NB) = AE3+AE1 
RHO I (NB ) = -AE2-AE4 
RHOD(NB) = RHOE(NB)+RHOI (NB) 

IF ( DV (MB ) .EO.O.O ) RHOD(NB) = 0.0 
RHOEOP ( NB ) = RHOE ( NB ) / SE 
RHODOP ( NB ) = RHOD ( NB ) /SE 
RHOIOP(NB) = RHO I ( NB ) / SE 
200 CONTINUE 

WNEPA = RHOEOP ( 1 ) 

COMPUTE EDV BY INTEGRATING RHOD (TRAPEZOIDAL RULE) 

202 CALL TRAP(RHOD»DV,NBB,EDV) 

DO 210 NB = 1 , NBB 
IF (EDV(NB).GT.O.O) EDV(NB)=0.0 
EDV(NB) = SORT ( C5*EDV ( NB ) ) 

210 CONTINUE 

EDVS = EOV(NBB) 

EDPZL = EDVS/ ( DVSRD/ AHDA ( K ) ) 

EE = EDVS 

IF (METS.EO.O) GO TO 400 

COMPUTE SC , EJS , TENT 

260 SC = SORT ( C4SEE ) 

PHD = -SC + PHI ( J) 

EJS(J) = 120.*TES*EXP (-PHD/AKTE ) 

TENT = .00 ISAM INK J 1 ( KOUNT ) , .12 ( KOUNT ) , J3 (KOUNT ) , J4( KOUNT ) , APJ ( J ) , 
1 AEJ ( J ) ) 

IF (KOUNT. EO.l) GO TO 340 
JA = J1 (KOUNT ) — J 1 (KUUNT-1 ) 

JB = J2(K0(JNT)— .J2IK0UNT-1) 

JC = J3( KOUNT) - J3 ( KOUNT- 1) 

CAT = AMAX1(JA,JB,JC) 

IF (CAT. LE. TENT ) GO TO 400 
GO TO 350 
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340 IF ( ABS ( EJ ( Jl-EJS ( J ) ) .LE. TENT ) GO TO 


400 


IF TEST IS NOT SATISFIED SET EJ = EJS 
RETURN TO COMPUTE A NEW DVS 

350 E J ( J ) = EJS(J) 

GO TO 38 

400 DO 405 NB = 1,N6B 

IF (EDV(NB).EQ.O.O) GO TO 403 
ESAVE(NB) = l.O/EDVINB) 

GO TO 405 

403 ESAVE(NB) = 0.0 
405 CONTINUE 

COMPUTE XDVS BY INTEGRATING 1.0/EDV (TRAPEZOIDAL RULE) 

CALL TRAP (ESAVE,DV,NBB,XSAVE) 

COMPUTE XL AM , TPN , CEN , TNE , PHAT , XAMTE , EAMPHZ 

XDVS = XSAVE(NBB) 

XLAM = XDVS/AMDA(K) 

TPN = 2.0*EPN(K) + APN 

CEN = ABS ( RHOI OP ( NBB ) ) + ABS ( RHOEOP ( NBB ) ) 

TNE = CEN + APN 

PHAT = (PHI (J)-PHIZ)/(AK*TE( I ) ) 

XAMTE = XDVS/ AMDATE ( K ) 

EAMPHZ = EDVS*AMDATE(K)/ (PHI (J ) -PH I Z ) 

COMPUTE XDVS 

DO 460 NB = 1 * NBB 
XDV(NB) = XDVS - XSAVE(NB) 

IF (NB.EO.NBB) XDV(NB) = 0.0 
460 CONTINUE 

IF ( IPRINT.EQ.l) GO TO 496 
IF ( METS.EQ.O ) GO TO 480 
WRITE (6,798) 

WRITE (6,800) AI,TE(I),PHI(J),EPN(K),TEP(L), 

1 , AMDATE ( K ) 

480 WRITE (6,810) 

WRITE (6,820) ( DV ( NB ) ,RHODOP ( NB ), RHOEOP ( NB ) , 

1 ( NB ) ,NB=1 ,NBB ) 

IF (METS.EO.l) GO TO 490 

WRITE (6,830) E J < J ) , EPJ ( J ) , P I J ( J ) *APJ (J),CJ(J),PP(J),EIJ(J),AEJ(J) 
1 ,SJA ( J ) ,SJ I ( J ) ,SJE< J ) ,SDJA( J) ,SOJE( J ) ,SDJ I ( J ) ,DVS( J ) ,XDVS ,APN 

2, XLAM,PHIZZ,EDVS 

3, TPM, CEN, TNE, PHAT, XAMTE, EAMPHZ ,WNEP A 

GO TO 496 

490 WRITE (6,831) E J ( J ) , EP J ( J ) » P I J < J ) » AP J ( J ) *C J ( J ) » PP ( J ) » E I J ( J ) ♦ AE J ( J ) 
1 ,SJA( J) ,SJI ( J ) ,SJE( J ) , SDJA ( J ) , SDJE ( J ) ,SDJ I ( J ) , DVS ( J ) , XDVS , APN 

2, XLAM, SC, PHIZ ,EDVS»DVSRD,DVSPZ,EDPZL»PHIZZ, DRDK 

3, TPN, CEN, TNE, PH AT, XAMTE, EAMPHZ ,WNEPA 
496 IF (METS.EO.l) GO TO 599 

C SAVE RICHARDSON - DUSHMAN VALUES OF DV , RHODOP , RHOEOP , EDV , 

C XDV , RHOIOP 

C 

DO 870 NB = 1 , NBB 
OVSAVE(NB) = DV(NB) 


TIP,AMDA(K),PTEST 
RHOI OP(NB),EDV(NB), XDV 
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RHDS(MB) = RHODUP ( MB ) 

RHES(NB) = RHOEQP ( NB ) 

EDS(NB) = EDV(NB) 

XDS(NB) = XDV(NB) 

RHIS(NB) = RHOIOP(NB) 

870 CONTINUE 
METS = 1 
IPRINT = 2 
GO TO 260 

CALL PLOTTING SUBROUTINE 

599 CALL PLOT 
970 CONTINUE 
980 CONTINUE 
990 CONTINUE 
99b CONTINUE 
1000 CONTINUE 
GO TO 1 

2 FORMAT ( 15/ ( 8E10.2 ) ) 

6 FORMAT ( 1H1, 1 OX , 90HP L ASMA TEST PR E S SURE ( P PT ) IS GREATER THAN THE V 
1AP0R PRESSURE OF THE PLASMA CHEM I C A L ( PT EST ) / 1 HO , 10X , 6HP PT = ,E12.5 
2,10X,8HPTEST = , E 1 2 . 5/ 1 HO, 1 OX , 5HT E = , F8 . 0 , 1 OX , 6HNE P = ,E8.1,10X,6 
3HTEP = ,F8.0) 

7 FORMAT ( 1H1 , 10X,99HDEBYE L ENGTH ( LAMBDA ) LONGER THAN MINIMUM MEAN F 
1REE PATH OF THE CHEM I CAL ( AMT ES T ) — COLLISIONAL C ASE / 1 HO , 10X , 10HL A 
2MBDA = , El 5,8, 10X»9HAMTEST = , E 1 5 . 8/ 1 HO , 10X , 5HTE = ,F8.0,6HNEP = , 
3E8.1 , 10X,6HTEP = ,F8.0) 

12 FORMAT (15) 

107 FORMAT ( 1 HO , 20X , 1 1 HDVS IS ZERO) 

109 FORMAT ( 1 HO , 20X , 20HDVS IS NEGATIVE STOP) 

798 FORMAT ( 1 H 1 , 5AX , 8HSC HOTTKY ) 

799 FORMAT ( 1 H 1 , ASX , 1 8HR I CHAR DSON-DUSHMAN ) 

800 FORMAT (1H0,2X,AHI = , F 5 . 3 , 5X , 5HTE = , F5 . 0 , 5X , 6HPH I = ,F5.3,5X,6HN 
1EP = , 1PE8.2,5X,6HTEP = , 1 PE R. 2 , 5 X , 6HT I P = , 0PF7. 1 , 5X , 9HL AMBDA = , 
21PE11.A/1H0,1X,5HPV = , 1 P E 1 3 . 6 , 5 X , 1 2HL AMBDA ( T E ) =,1PE11.A> 

801 FORMAT (1H1,2X,AHI = , F 5 . 3 , 5 X , 5HT E = , F5 . 0 , 5 X , 6H PH I = ,F5.3,5X,6HN 
1EP = ,1PE8.2,5X,6HTEP = , 1 P E 8 . 2 , 5X , 6HT I P = , 0 PF 7 . 1 , 5X , 9HL AMBDA = , 
21PE11.A/1H0,1X,5HPV = , 1PE13.6,5X, 12HLAMBDA( TE ) =,1PE11.A) 

810 FORMAT ( 1 HL , 1 OX , 2HD V , 12 X , 6HND ( D V ) , 1 2 X , 6HNE ( DV ) , 1 2 X , 6HN I ( DV ) , 1 2 X , 5H 
2 E ( DV ) , 12X,5HX(DV)/1H0) 

820 FORMAT ( OPF 16 . 5 , 1 PE 1 9 . 6 , 1 P2E 1 8 . 6 , 1 P E 1 7 .6 , 1 PE 1 8 . 6 ) 

830 FORMAT ( 1 HL , 2 X , 9H JEE = , 1 PE 1 3 . 6 , AX , 9HJ E P = , 1 P E 1 3. 6 , AX , 6HJ I P 

1 = , 1PE13.6,AX,9HJAP = ,1PE13.6/ 

23X ,9HJ = ,1PE13.6,AX,9HPP = , 1 P E 1 3. 6 , AX , 6H J I E = ,1PE13.6, 

3AX.9HJAE = , 1PE13.6/ 

A3X , 9H J A = ,1PE13.6,AX,9HJI = , 1 P E 1 3. 6 , AX , 6HJ E = ,1PE13.6, 

5AX,9HJA/JAP = ,1PE13.6/ 

63X,9HJE/JEP = , 1PE13.6,AX,9HJ 1/ J I P = , 1 P E 1 3. 6 , AX , 6HDVS = ,0PF8.5,9 

7 X , 9HXDV S = ,1PE13.6/3X,9HNAP = , 1 PE 1 3 . 6 , AX , 9HXD/ L AM = ,1PE13. 

86 , AX, 6HPHZ Z= , 0PF8 . 5 , 9 X , 9HE DVS = ,1PE13.6/ 

9 3X , 9HNT P = , 1PE13.6,AX,9HNCE = , 1 P E 1 3 .6 , AX , 6HNTE = ,1PE13.6, 

1 AX , 9HRD/KTE = , 1 PE 13, 6/3X , 9HX/LMTE = , 1 P E 1 3. 6 , AX , 9HE L T/ RD = ,1PE13 

2.6,AX,6HNEPA= ,1PE13.6) 

831 FORMAT ( 1 HL , 2 X , 9HJEE = , 1 PE 1 3 . 6 , AX , 9H J E P = , 1 P E 1 3. 6 , AX , 6HJ I P 

1 = ,1PE13.6,AX,9HJAP = ,1PE13.6/ 

2 3X , 9HJ = ,1PE13.6,AX,9HPP = , 1 PE 1 3 . 6 , AX , 6HJ I E = ,1PE13.6, 

3 AX , 9HJ A E = , 1PE13.6/ 

A3X , 9HJ A = , 1PE13.6, AX, 9HJ I = , 1 P E 1 3. 6 , AX , 6HJ E = ,1PE13.6, 

5AX »9HJ A/ JAP = ,1PE13,6/ 
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63X,9HJE/JEP = ,1PE13.6,4X,9HJ I/JIP = , 1 PE 1 3. 6 , 4X , 6HDVS = ,OPF8.5,9 
7 X , 9HXDV S = ,1PE13.6/3X,9HNAP = , 1 PE 1 3 .6 , 4X , 9HXD/ L AM = ,1PE13. 

86,4X,6HSC = ,1PE13.6,4X,9HPHZ = ,0PF8.5/ 

93X , 9HEDVS = , 1 PE 1 3 . 6 , 4X , 9HDVSRD = , 1 P E 13 . 6 , 2 7X , 9HDV S /RD = ,1PE1 
13.6/ 

23X * 9HELM/ RD = , 1 PE 1 3 . 6 , 4X , 9HPHZ Z = , 1 PE 13 . 6 , 2 7 X r 9HDRD/KT = ,1PE1 

33.6/ 

43X » 9HNTP = , 1PE13.6,4X,9HNCE = , 1 P El 3 . 6 , 4X , 6HNT E = ,1PE13.6, 

54X , 9HRD/KTE = , 1 PE 13. 6/3X , 9HX/LMTE = , 1 PE 1 3 . 6 , 4X , 9HEL T/ RD = ,1PE13 
6.6,4X,6HNEPA= ,1PE13.6> 

840 FORMAT ( 1 HL , 2 X , 9HRH0 ( DV ) = ,1PE13.6,3H + , 1 PE 1 3 . 6 , 6H*DV + ,1PE13.6, 
19H=!=DV=?*2 + ,1PE13.6,9H*DV**3 + , 1 PE 1 3 .6 1 6H*DV**4 ) 

850 FORMAT (1H1) 

860 FORMAT ( 1HL/ 1 HL , 2 X , 4H I = , F5 . 3 , 5X , 5HTE = , F5 .0 , 5 X T 6HPH I = ,F5.3,5X 
1 » 6HNEP = ,1PE8.2,5X,6HTEP = , 1 P E8 . 2 , 5X , 6HT I P = , 0PF7. 1 , 5X, 9HL AMBDA 
2 = ,1PE11.4/1H0,1X,5HPV = , 1 P E 1 3 .6 , 5X , 1 2HL AMBDA ( T E ) =,1PE11.4) 

END 

$ I BFTC PLOTA 

SUBROUTINE PLOT 

COMMON /HP/ XDV ( 50 ) , XDSI50) ,RHODOP( 50 ) ,RHDS( 50 ) , EDV ( 50 ) , E DS ( 50 ) t 
1DV( 50) , DVSAVEt 50) ,RHOEOP( 50) , RHES ( 50 ) , RHO I OP ( 50 ) , RH I S ( 50) ,NBB 
2, I WRITE 

DIMENSION KKK ( 14 ) , P ( 10 ) , Z ( 100 ) , Z A ( 1 00 ) , ZB ( 1 00 ) , ZC ( 1 00 ) , Z D ( 100 ) , 

1 Z E ( 100 ) 

599 ND = 0 


DO 600 NB 

= 1 1 NBB 

NE = NBB 


- ND 

NET = 2* 

NBB - ND 

ND = ND 

+ 

1 

Z { NE ) = 


XDV(NB) 

Z ( NET ) = 


XDS(NB) 

Z A ( NE ) 

= 

RHODOP i NB ) 

ZA(NET) 

= 

RHDS (NB ) 

Z B { NE ) 

= 

EDV(NB) 

ZB(NET) 

= 

EDS( NB) 

ZC l NE ) 

= 

DV(NB) 

ZC ( NET ) 

= 

DVSAVE(NB) 

Z D { NE ) 

= 

RHOEOP ( NB ) 

ZD(NET) 

- 

RHES(NB) 

ZE(NE) = 


RHOIOP(NB) 

ZE(NET) 

= 

RH I S < NB ) 


600 CONTINUE 

PI 1 ) = 5.0 
KKK ( 1 ) = 64 

KKK I 2 ) = 2 
KKK I 3 ) = NBB 
KKK ( 5 ) = NBB 
NB2 = 2*NBB 

CALL SCALE ( NB2 , Z A , KR STR ) 
CALL PLOT MY <Z,ZA,KKK,P) 
WRITE (6,602) KRSTR 
ND = 0 

DO 607 NB = 1 , NBB 
NE = NBB - NO 
NET = 2*NBB - ND 
Z(NE) = XDV(NB) 

Z ( NET ) = XOS(NB) 

NO = NO + 1 
607 CONTINUE 

CALL SCALE ( NB2 , Z D , KRSTR ) 



CALL PLOTMY (Z,ZD,KKK,P) 
WRITE (6,603) KRSTR 
ND = 0 

DO 60 L NB = 1 , NBB 
NE = NBB - NO 
NET = 2*N8B - ND 
ND = ND + 1 
Z ( NE ) = XDV(NB) 

Z ( NET ) = XDS(NB) 

601 CONTINUE 

CALL SCALE ( NB2 , Z E , KRSTR ) 
CALL PLOTMY(Z,ZE,KKK,P) 
WRITE (6,610) KRSTR 
ND = 0 

DO 611 NB = 1 , NBB 
NE = NBB-ND 
NET = 2*NBB-ND 
ND = ND+1 
Z ( NE ) = XDV(NB) 

Z ( NET ) = XDS(NB) 

611 CONTINUE 

CALL SCALE ( NB2 , ZB , KRSTR ) 
CALL PLOTMY <Z,ZB,KKK,P) 
WRITE (6,604) KRSTR 
ND = 0 

DO 606 NB = 1 »NBB 
NE n NBB - ND 
NET = 2*NBB - ND 
ND = ND + 1 
Z ( NE ) = XDV(NB) 

Z ( NET ) = XDS(NB) 

606 CONTINUE 

CALL SCALE ( NB2 , ZC , KRSTR ) 
CALL PLOTMY (Z,ZC,KKK,P) 
WRITE (6,605) KRSTR 
IF ( IWRITE.EO.l ) RETURN 
P ( 1 ) = NBB 
KODE = 64 
ND = 0 

DO 500 NB = 1 , NBB 
NE = NBB - ND 
Z ( NE ) = XDV(NB) 

ZA(NE) = RHODOP ( NB ) 

ZD(NE) = RHOEOP ( NB ) 

ZB(NE) = EDV(NB) 

ZC(NE) = DV(NB) 

ZE(NE) = RHOIOP(NB) 

ND = NO + 1 
500 CONTINUE 

WRITE (6,501) 

CALL SCALE ( NBB , Z A , KRSTR ) 
CALL PLOTXY ( Z , Z A ,KODE , P ) 
WRITE (6,502) KRSTR 
WRITE (6,501) 

CALL SCALE < NBB , ZD , KRSTR ) 
CALL PLOTXY ( Z , ZD , KODE , P ) 
WRITE (6,505) KRSTR 
WRITE (6,501) 

CALL SCALE ( NBB , Z E , KR STR ) 
CALL PLOTXY ( Z , Z E , KODE , P ) 
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WRITE (6,510) KRSTR 
WRITE (6,501) 

CALL SCALE ( NBB , ZB , KR STR ) 

CALL PLOTXY ( Z , ZB , KOOE , P ) 

WRITE (6,503) KRSTR 
WRITE (6,501) 

CALL SCALE ( NBB , ZC , KRSTR ) 

CALL PLOTXY ( Z , ZC ,KOI)E , P ) 

WRITE ( 6 , 504 ) KRSTR 

501 FORMAT ( 2HPT ) 

502 FORMAT ( 2HPL , 47 X , 8HN0 ( X 10** , I 3 , 6H ) VS X) 

503 FORMAT ( 2HPL , 47X , 7HE ( X 10** , 1 3 , 6H ) VS X) 

504 FORMAT ( 2HPL , 47 X , 8HDV ( X 1 0** , 1 3 , 6H ) VS X) 

505 FORMAT ( 2HPL ,47X , 8HNE ( X 10** , I 3 , 6H ) VS X) 

510 FORMAT < 2HPL , 47 X , 8HN I ( X 1 0** , I 3 , 6H ) VS X) 

602 FORMAT ( 2HPL , 47X , 8HN0 ( X 10** , I 3 , 6H ) VS X/ 

12HPL,44X,20H+ R I CHARDS0N-DUSHMAN/2HPL , 44X , 1 OH* SCHOTTKY) 

603 FORMAT ( 2HPL ,47X , 8HNE ( X 10** , I 3 , 6H ) VS X/ 

12HPL,44X, 20H+ R I CHARDSON-DUSHMAN/ 2HPL , 44X , 1 OH* SCHOTTKY) 

604 FORMAT ( 2H PL , 47 X , 7HE ( X 10** , 1 3 , 6H ) VS X/ 

12HPL,44X,20H+ R I CHARDSON-DUSHMAN/ 2HPL , 44X , 1 OH* SCHOTTKY) 

605 FORMAT ( 2HPL ,47X , 8HDV ( X 10** , I 3 , 6H ) VS X/ 

12HPL,44X,20H+ R ICHARDSON-DUSHMAM/ 2HPL , 44X , 10H* SCHOTTKY) 

610 FORMAT ( 2HPL ,47X , 8HN I ( X 10** , I 3 , 6H ) VS X/ 

12HPL,44X,20H+ R I CHARDSON-DUSHMAN/ 2HPL , 44X, 10H* SCHOTTKY) 
RETURN 
END 

$ I BFTC VINE 

SUBROUTINE VIN ( NV , DVV , T , ANS ) 

COMMON /MV/ PI, EM, AK, AM 

IF NV = 1 COMPUTATIONS FOR ELECTRONS 

IF NV = 2 COMPUTATIONS FOR ATOMS AND IONS 

THE PROPER T AND DV IS SPECIFIED IN MAIN PROGRAM 

IF (NV.EO.l) GO TO 10 
C = SORT (2.0*AK*T/(PI*AM) ) 

60 TO 20 

10 C = SORT (2.0*AK*T/(PI*EM) ) 

20 XE = SORT (DVV/(AK*T) ) 

Y = ERF(XE) 

ERFC = 1.0-Y 

ANS = C*EXP(-DVV/(AK*T) J/ERFC 

RETURN 

END 

I BFTC TRAPE 

SUBROUTINE TRAP ( X , DV , NBB , ANS ) 

DIMENSION X(40) »DV( 50) , ANSI 40) 

SUBROUTINE TO INTEGRATE TR AP EZO I DAL L Y 

H = DV ( 2 ) /2.0 

SUME = 0.0 

DO 60 NB = 1 ,NBB 

IF (NB.EO.l ) GO TO 50 

SUME = SUME+X(NB-1 )+X(NB) 

ANS(NB) = SUME*H 
GO TO 60 

50 ANS ( N8 ) = 0.0 
60 CONTINUE 



$DATA 


RETURN 

END 


1 

2.4E+3 

1 

4.0E+0 

1 

l.E+13 

1 

2.5E+3 

1 

2.4E+3 
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APPENDIX B 


ELECTRON SHEATH PROGRAM LISTING 


$ I BFTC (ELECT 

A PLANAR ELECTRON SHEATH BETWEEN AN EMITTER AND A 
NEAR-EQUILIBRIUM PLASMA 

COMMON /MV/ PI, EM, AK, AM 

COMMON /MP/ XDV( 50) ,XDS< 50) ,RHGDUP<50) ,RHDS(50) , EDV ( 50 ) , EDS ( 50 ) , 
1DV( 50 ) ,DVSAVE( 50) ,RHOEOP( 50 ) , RHES ( 50 ) , RHO I OP ( 50 ) ,RHIS ( 50) ,N8B, 
2IWRITE 

DIMENSION TE< 10) ,PHI ( 10) ,EPN( 10) ,TEP( 10) ,TIPP( 10) ,EJ ( 10) ,EPJ( 10) , 
1PIJ (10) ,APJ ( 10) ,PP< 10) ,OVS( 10) ,EI J ( 10) , AEJ< 10) ,SJE( 10) ,SJI ( 10) , 
2SJA( 10) ,CJ ( 10) ,SDJA( 10) ,SOJI ( 10 ) , SDJ E ( 1 0 ) , RHOD ( 50 ) ,RH01 ( 50) , 

3RH06 ( 50 ) , ES AVE ( 50 ) ,XSAVE(50) , E I J S ( 10 ) , AMDA ( 1 0 ) , 

4AMDATE ( 10 ) 

A I , IONIZATION POTENTIAL 

301 A I = 3.893 
AK = 8.617E-5 
AKT = 2.0*AK 

EM = .511E+6/( (2.998E+10)**2 ) 

PI = 3.1^159 
PI EM = PI*EM 
SE = 1.602E-19 

MA , MASS OF ATOM 

302 AM = 931.478E+6*132.9/( (2.998E+10>**2) 

P I AM = P I *AM 
API = 6.6256E-34/SE 
A P2 = AM 

Cl = SQRT(AK/(4.0*PI*.511E+6*2.82E-13)) 

C2 = AP 1 #$2 

C3 = 82. 06*760. /6.023E+23 
C5 = -72.0*PI*1.0E+11 

IWRITE = 1 FOR SCHOTTKY OUTPUT 

IWRITE = 0 FOR RICHARUSON-DUSHMAN AND SCHOTTKY OUTPUT 
READ (5,12) IWRITE 

READ INPUT VALUES OF TE , PHI , EPN , TEP , TIP 

1 READ (5,2) I I , (TE( I ) , 1 = 1, II ) 

READ (5,2) JJ,(PHI(J),J=1,JJ) 

READ (5,2) KK, ( EPN(K) ,K=1 ,KK) 

READ (5,2) LL, (TEP(L) » L = 1 » LL ) 

READ (5,2) LA, (TIPP(LB) ,LB=1»LA) 

DO 1000 I = 1,11 

COMPUTE PTEST , VAPOR PRESSURE 
C 

303 PTEST = 10.**(-3920.38/TE( I ) -. 5 1 97 8 1*ALGG 10 < TE < I ) )+10. 71914)/ 
1133.322 

AKTE = AK*IE(I) 

TAKTE = 2 .0* AKT E 
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TES = TE( I ) *TE ( I ) 

SB = SORT ( TAKTE/P I EM) 

DO 960 LB = 1,LA 
TIP = TIPP(LB) 

AKT P = AK*TIP 
TAKTIP = AKT*T I P 

SET TAP = TIP FOR THESE SOLUTIONS 

TAP = TIP 

TAKTAP = AKT*TAP 

SI = SORT ( TAKTIP/ PI AH) 

S 1 A = SORT (TAKTAP/PIAM) 

DO 990 L = 1 , L L 
AKTEP = AK*TEP(L) 

TAKT6P = 2.0*AKTEP 
SA = SORT ( TAKTEP/P I EM ) 

CL1 = TAKTEP*PIEM 
DO 980 K = 1 * KK 
SEE = SE*EPN(K) 

COMPUTE PHIZZ , AM DA , AMDATE , APN , PPT , AMTEST 

PHIZZ = -AKTE*ALDG( (SEE/(240.*TES) >*SB) 

AMD A ( K ) = C1*S0RT(TEP(L)/EPN(K> ) 

AMDATE (K) = C1*SQRT ( TE ( I > / EPN ( K ) > 

APN = ( EPN(K)*C2 )*( EPN(K)*( API ) )*EXP( A I / AKTEP ) / ( CL 1** 1 . 5 ) 

PPT = C3*(2.0*EPN(K)*TIP+APN*TAP) 

AMTEST, MINIMUM MEAN FREE PATH 

304 AMTEST = 1.0E+12/APN 

TEST ON AMDA AND AMTEST 

IF ( AMDA(K ) .LT. AMTEST ) GO TO 9 

WRITE (6,7) AMDA(K) , AMTEST, TE< I ), EPN(K ), TEP< L ) 

GO TO 980 

TEST ON PPT AND PTEST 
9 IF (PPT -PTEST) 40,40,6 

6 WRITE (6,6) PPT , PTEST ,TE( I ) ,EPN(K) »TEP(L) 

GO TO 980 
40 CONTINUE 

DO 970 J = 1,JJ 
METS = 0 
IPRINT = IWRITE 
3 PHD = PHI ( J ) 

COMPUTE EPJ , PIJ , A P J , PP 

EPJ(J) = SEE*.S*SA 

PIJ(J) = SEE*.5*SI 

APJ(J) = SE*APN*SI A/2.0 

PP(J) = C3*(EPN(K)*TEP(L)+EPN(K)*TIP+TAP*APN> 

IF ( IWRITE. EO.l. OR .METS.EO.l) GO TO 38 
WRITE (6,799) 

8 WRITE (6,800) A I , T E ( I ) , PH I ( J ) , EPN ( K ) , T EP ( L ) , T I P , AMDA ( K ) , PTEST 

1, AMDATE (K) 
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38 CONTINUE 
COMPUTE EJ 

E J ( J ) = 120.*TES*EXP( (-PHD ) /AKTE ) 

COMPUTE DVS 

OVS(J) = AKTE*AL0G(SEE*SB/(2.0*EJ( J) ) ) 

IF ( METS • EQ .0 ) DVSRD = DVS(J) 

PHIZ = PHI ( J ) - DVSRD 
DRDK = ABS(DVSRD/AKTEP) 

DVSPZ = DVS ( J )/ DVSRD 
115 IF ( DVS l J ) >102,106,108 
108 WRITE (6,109) 

GO TO 970 
106 WRITE (6,107) 

GO TO 970 

COMPUTE E I J , AEJ , SJE , SJI , SJA , CJ , SDJI , SDJE , SDJA 

102 EIJ(J) = (APJ(J)+PIJ(J)*EXP(DVS( J)/AKTP ) ) / < 1 .0+2 .0*EXP ( ( 

1 A I-PHD ) / AKTE ) ) 

AEJ(J) = (APJ (J)+PIJ(J)*EXP(DVS(J)/AKTP ) ) / ( 1 .0+ . 5* EX P ( - ( 

1 ( A I-PHD ) /AKTE ) ) ) 

SJE(J) = EJ(J)*EXP(DVS(J)/AKTE >-EPJ(J> 

S J I ( J ) = P I J ( J)*EXP(DVS( J)/AKTP ) - EIJ(J) 

SJA(J) = AEJ(J) - APJ(J) 

C J ( J ) = SJE(J) + SJI (J) 

SDJI(J) = S J I ( J ) / PIJ(J) 

SDJE(J) = SJE(J) / EPJ(J) 

SDJA(J) = SJA(J) / APJ(J) 

IN = 20 

COMPUTE - IN - VALUES OF DV FROM 0.0 TO DVS 
IN = 20 

112 DVI = DVS ( J ) /FLOAT ( IN) 

DV(1) = 0.0 
DO 120 NB = 1 , 1 N 
NBB = NB + 1 
DV(NBB) = DV(NB) + DVI 
120 CONTINUE 

COMPUTE RHOEOP , RHOIOP , RHODUP 
ERF IS THE ERROR FUNCTION SUBROUTINE 
VIN IS A SUBROUTINE TO CALCULATE SORT ( 2KT/P I*M ) *EXP (-DV/KT ) / 
( 1.0-ERF(S0RT(DV/KT) ) ) 


121 DO 200 NB = 1 , NBB 

X = SORT (-DV(NB)/AKTE) 

ERA = ERF ( X ) 

X = SORT(-(DVS(J)“DV(NB) )/AKTP) 
ERB = ERF ( X ) 

DV P = DV(NB) 

IF (DVP.EQ.O.O) GO TO 160 
CALL VIN ( 1,DVP,TEP(L) ,ANS) 

AVI = ANS 
GO TO 161 

160 AVI = SORT (TAKTEP/PIEM) 
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161 CONTINUE 

DVP = UVStJ) - OV(NH) 

IE (DVP. EQ. 0.0) GU TO 165 
CALL VIN ( 2 » DVP »TE ( I ) ,ANS) 

A V2 = ANS 
GO TO 166 

165 AV2 = SORT (TAKTE/PIAM) 

166 CONTINUE 

EA = EXP((DVS(J)-DV(NB))/AKTE) 

AE1 = E J ( J ) *FA* ( 1 • 0+ERA ) / SB 
AE2 = EPJ(J)/AV1 
AE3 = E I J ( J ) / AV2 
FB = EXP(DV(NB)/AKTP) 

AF4 = PIJ(J)*FB*(1.0+I:RB)/SI 
RHOE(NB) = AE1+AE2 
RHOI(NB) = -AE3-AE4 
RHOD(NH) = RHOE(NB)+RHUHIMB) 

IE (DV(NB).EQ.O.O) RHOD(NB) = 0.0 
RHOEOP(NB) = RHUE ( MB ) / SE 
RHODOP(NB) = RHDD(NB)/SE 
RHGIOP(NB) = R HO I (NB)/SE 
200 CONTINUE 

WNIPA = RHOIOP(l) 

COMPUTE EDV BY INTEGRATING RHUL) (TRAPEZOIDAL RULE) 

202 CALL TRAP ( RHOOtDVt NBBt EDV ) 

DO 210 NB = IfNBB 

IE (EDV(NB).GT.O.O) EDV(NB) = 0.0 
EDV(NB) = -S0RT(C5*E!>V(NB) ) 

210 CONTINUE 

EDV S = EDV(NBB) 

EDPZL = EDVS/ ( DVSRD/ AMOA ( K ) ) 

EE = EDVS 

IF (METS.EO.O) GO TO 402 
COMPUTE SC t E I JS 

260 SC = S0RT(-.511E+6*2.82E-13*EE) 

PHD = SC + PHI ( J ) 

EIJS(J) = (APJ(J)+PIJ (J)-EXP(DVS(J J/AKTP ) ) / ( 1 .0+2 . 0*EXP ( < A I 

1PHD) / AKTE ) ) 

TENT = • 001 -AM INI ( E J ( J ) , EP J ( J ) , P I J ( J ) t E I J ( J ) , A PJ ( J ) , A E J ( J ) ) 

IE (ABS(EIJ(J ) - E I JS ( J ) ) . L E . TENT ) GO TO 402 

IF TEST IS NOT SATISFIED SET EIJ = EIJS AND RETURN 
TO COMPUTE A NEW DVS 

E I J ( J ) = E I JS ( J ) 

GO TO 38 

402 DO 405 m = 1 j NBB 

IF (EI)V(NB) .EQ.O.O) GO TO 403 
ESAVE(NB) = 1 • 0/ EDV ( NB ) 

GO TO 405 

403 ESAVE(NB) = 0.0 
405 CONTINUE 

COMPUTE XDVS BY INTEGRATING 1.0/EUV (TRAPEZOIDAL RULE) 

CALL TRAP ( ESAVEtDVtNBBtXSAVE ) 
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XDVS = XSAVE ( NBB ) 


C 

C COMPUTE XL AM , TPN , CEN , TNE , PHAT , XAMTE , EAMPHZ 

C 

XL AM = XDVS/ AMOA ( K ) 

TPN = APN + 2 .0*EPN ( K ) 

CEN = RHOEOP(NBB)+ABS(RHOIOP(NBB) ) 

TNE = CEN + APN 

PHAT = (PHI (J)-PHIZ)/ ( AK*TE ( I ) ) 

XAMTE = XDVS/ AMDATE ( K ) 

EAMPHZ = EDVS*AMDATE(K)/(PHI (J)-PHIZ) 

C 

C COMPUTE XDV 

C 

DO 460 NB = 1 » NBB 
XDV(NB) = XDVS - XSAVE ( NB ) 

IF (NB.EO.NBB) XDV(NB) = 0.0 
460 CONTINUE 

IF ( IPRINT.EO.l ) GO TO 496 
IF (METS.EO.O) GO TO 480 
WRITE ( 6 » 798 ) 

WRITE (6,800) AI ,TE( I ) ,PHI (J) ,EPN(K) ,TEP(L ) ,T IP, AMDA(K) ,PTEST 
1 , AMDATE ( K ) 

480 WRITE (6,810) 

WRITE (6,820) ( DV ( NB ) , RHODOP ( NB ) , RHUEOP ( NB ) , RHO 1 OP ( NB ) , EDV ( NB ) 

1 , XDV ( NB ) ,NB=1,NBB) 

IF (METS.EQ.l ) GO TO 490 

WRITE (6,830) E J ( J ) , FPJ ( J ) , P I J ( J ) , AP J ( J ) ,C J ( J ) , PP ( J ) , E I J ( J ) , AE J ( J ) 
1 , SJ A ( J ) ,SJI (J) ,SJE(J) , SDJ A ( J ) ,SDJE(J ) ,SDJI (J ) ,DVS(J) ,XDVS ,APN 

2, XLAM,PHIZZ,EDVS 

3, TPN, CEN, TNE, PHAT, XAMTE, EAMPHZ, WN I PA 
GO TO 496 

490 WRITE (6,831) E J ( J ) , EP J ( J ) , P 1 0 ( J ) , AP J ( J ) ,C J ( J ) , PP ( J ) , E I J ( J ) , AE J ( J ) 
1,SJA(J),SJI(J),SJE(J)»SDJA(J),SDJE(J), SDJ I ( J ) , DVS ( J ) , XDVS ,APN, 

2 XL AM, SC, PHIZ , EDVS , DVSRD, DVSPZ , EDPZL , PH I Z Z , DRDK 
3, TPN, CEN, TNE, PHAT , XAMTE , EAMPHZ ,WN IP A 
496 IF (METS.EO.l) GO TO 599 

SAVE RICHARDSON - DUSHMAN VALUES FOR OV , RHODOP , RHOEOP , 

RHQIOP , EDV , XDV 

DO 870 NB = 1 , NBB 
DVSAVE ( NB ) = DV(NB) 

RHDS(NB) = RHODOP ( NB ) 

RHES(NB) = RHOEOP ( NB ) 

EDS(NB) = EDV(NB) 

XDS(NB) = XDV(NB) 

RHIS(NB) = RHOIOP(NB) 

870 CONTINUE 
METS = 1 
IPRINT = 2 
GO TO 260 

CALL PLOTTING SUBROUTINE 

599 CALL PLOT 
970 CONTINUE 
980 CONTINUE 
990 CONTINUE 
960 CONTINUE 
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1000 CONTINUE 
GO TO 1 

2 FORMAT ( I5/(8E10.2) ) 

6 FORMAT (1H1»10X»90HPLASMA TEST PRESSURE ( PPT ) IS GREATER THAN THE V 
1AP0R PRESSURE OF THE PLASMA CHEM IC AL ( PTEST ) / 1H0, 10X , 6HPPT = ,E12.5 
2,10X,8HPTEST = , E 1 2 . 5 / 1H0 , 10X , 5HTE = , F8 .0 , 1 OX , 6HNE P = ,E8.1,10X,6 
3HTEP = » F8. 0 ) 

7 FORMAT ( 1H1 , 10X .99HDEBYE LENGTH ( LAMBDA ) LONGER THAN MINIMUM MEAN F 
1REE PATH OF THE CHEM I C AL I AMTEST ) — COLLISIONAL CASE/ 1H0, 10X , 10HLA 
2MBDA = , E15.8 r 10X,9HAMTEST = » E 1 5 . 8 / 1H0 , 10X , 5HTE = ,F8.0,6HNEP = , 
3E8.1, 10X,6HTEP = ,F8.0) 

12 FORMAT (15) 

107 FORMAT ( 1H0 , 20X , 1 1 HDV S IS ZERO) 

109 FORMAT ( 1H1 ,20X ,20H0VS IS POSITIVE STOP) 

798 FORMAT < 1H1 , 54X , 8HSCH0TTKY ) 

799 FORMAT <1H1»48X»18HRI CHAROSON— OUSHMAN ) 

800 FORMAT <lH0,2X f 4HI = , F5 . 3 , 5X , 5HTE = , F5 .0, 5X, 6HPHI = ,F5.3,5X,6HN 
1EP = , 1PE8.2 ,5X,6HTEP = , 0PF7 . 1 , 5X , 6HT I P = , 0PF7 . 1 , 5 X , 9HLAMB0A = , 
21PE11.4/1H0,1X,5HPV = , 1 PE 1 3 . 6 , 5X , 1 2HL AMBDA ( TE ) =,1PE11.4) 

801 FORMAT <1H1,2X,4HI = , F5. 3, 5X ,5HTE = , F 5. 0 , 5X , 6HPH I = ,F5.3 f 5X,6HN 
1EP = ,1PE8.2,5X,6HTEP = , OPF 7 . 1 , 5 X , 6HT I P = , 0PF7 . 1 , 5X , 9HL AMBDA = , 
21PE11.4/1H0,1X,5HPV = ,LPE13.6,5X,12HLAMBDA(TE) =,1PE11.4) 

810 FORMAT ( 1HL , 10X , 2HDV , 1 2 X , 6HND ( DV ) , 1 2X , 6HNE ( DV ) , 1 2 X , 6HN I ( OV ) , 12X,5H 
1E(DV),12X,5HX(DV)/1H0) 

820 FORMAT ( OPF 1 6 . 5 , 1 PE 19 . 6 , 1 P2E 1 8 . 6 , 1 PE 1 7 . 6 , 1 PE 1 8 . 6 ) 

830 FORMAT ( 1 HL , 2 X , 9HJ E E = , 1PE1 3 .6 ,4X ,9HJEP = , 1 PE 1 3 .6 , AX , 6H J I P 

1 = ,1PE13.6, 4X.9HJAP = ,1PE13.6/ 

23X T 9HJ = , 1PE13.6,4X,9HPP = , 1 PE 13 .6 , AX , 6H J I E = .1PE13.6, 

34X,9HJAE = * 1 PE 13 .6/ 

43X,9HJA = , 1PE13.6,4X,9HJ I = , 1 PE 1 3 . 6 , AX , 6H J E = .1PE13.6, 

54X * 9HJ A / J AP = t 1PE13.6/ 

63X » 9HJE/ JE P = ,1PE13.6,4X,9HJI/JIP = , 1 PE 1 3 . 6 , 4X * 6HDV S = f 0PF8.5»9 
7X.9HXDVS = ,1PE13.6/3X,9HNAP = , 1 PE 1 3. 6 , 4X , 9HXD/L AM = ,1PE13. 

86,4X,6HPHZZ= ,0PF8 . 5 , 9X , 9HEDV S = ,1PE13.6/ 

93X,9HNTP = , 1PE13.6,4X,9HNCE = , 1PE 1 3 . 6 , 4X , 6HNTE = ,1PE13.6, 

14X.9HR0/KTE = , 1 PE 1 3 . 6/ 3X ,9HX/ LMTE = , 1 PE 1 3 . 6 , 4X , 9HEL T /RD = ,1PE13 

2.6,4X,6HNIPA= .1PE13.6) 

831 FORMAT ( 1HL , 2 X , 9HJ EE = , 1PE1 3 .6 ,4X r 9HJEP = , 1PE 1 3 .6 , 4X , 6H J I P 

1 = , 1PE13.6,4X,9HJAP = ,1PE13.6/ 

23X,9HJ = t lPE13.6,4X,9HPP = , l PE 1 3 .6 , 4X , 6H J I E = ,1PE13.6, 

34X.9HJAE = * 1 PE 13 .6/ 

43X.9HJA = , lPE13.6»4X f 9HJ I = , 1 PE 13 .6 , 4X , 6H JE = ,1PE13.6, 

54X,9HJA/JAP = » 1 PE 1 3 . 6/ 

63X * 9HJE / J EP = ,1PE13.6,4X,9HJI/JIP = r 1 PE 1 3 . 6 , 4X , 6HDVS = ,0PF8.5,9 

7 X » 9HXDV S = ,1PE13.6/3X,9HNAP = , 1 PE 1 3 . 6 , 4X , 9HXD/L AM = ,1PE13. 

86 * 4X 1 6HSC = ,1PE13.6,4X,9HPHZ = ,0PF8.5/ 

9 3X * 9HEDV S = , 1 PE 13 . 6 , 4X , 9HDVSRD = t 1PE 1 3. 6 , 2 7X, 9HDV S/ RD = ,1PE1 

13.6/ 

23X.9HELM/RD = , 1 PE 1 3 . 6 , 4X , 9HPHZZ = , 1 PE 1 3. 6 , 27 X t 9HDRD/K T = ,1PE1 

33.6/ 

43X * 9HNTP = , 1PE13.6,4X,9HNCE = , 1 PE 1 3 . 6 , 4X , 6HNTE = ,1PE13.6, 

54X,9HRD/KTE = , 1 PE 1 3 . 6/ 3X , 9HX/ LMT E = , 1 PEI 3 . 6 , 4X , 9HELT /RD = ,1PE13 

6.6,4X,6HNIPA= ,1PE13.6» 

840 FORMAT ( 1 HL , 2 X , 9HRH0 < DV ) = ,1PE13.6,3H + , 1 PE 13 . 6 , 6H*DV + *1PE13.6, 
19H*DV=x*2 + , 1PE13.6,9H*DV**3 + , 1 PE 1 3. 6 , 6H*DV**4 ) 

850 FORMAT (1H1) 

860 FORMAT ( 1 HL/ 1HL , 2 X , 4H I = , F5 . 3 , 5X , 5HTE = , F5. 0 , 5X , 6HPH I = t F5.3,5X 
1 r 6HNE P = ,1PE8.2,5X,6HTEP = , 0PF7 . 1 , 5 X , 6HT I P = , 0PF7 . 1 , 5X , 9HLAMBDA 

2 = ,1PE11.4/1H0,1X,5HPV = , 1 PE 1 3 . 6 , 5X , 12HL AMBDA ( TE ) =,1PE11.4> 

END 
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$ I B P TC PLOTA 

SUBROUTINE PLOT 

COMMON /MP/ XDV< 50) ,XDS< 50) ,RHODOP( 50) ,RHDS( 50) ,EDV(50) ,EDS( 50) , 
1DV( 50) ,DVSAVE(50) ,RH0E0P(50) , RHES ( 50 ) ,RHOIOP ( 50 ) ,RH I S ( 50 ) ,NBB, 

21 WRITE 

DIMENSION KKM 14) f P( 10) , Z < 100 ) , Z A ( 1 00 ) , Z B ( 100 ) , ZC ( 1 00 ) , Z D ( 1 00 ) , 

1 Z E { 100) 

599 NO = 0 

DO 600 NB = 1 , NBB 
NE = NBB - NO 
NET = 2*NBB - ND 
ND = ND + 1 
Z ( NE ) = XDV(NB) 

Z(NET) = XDStNB) 

ZA(NE) = RHODOP(NB) 

Z A ( NET ) = RHDS(NB) 

Z B ( NE ) = EDV(NB) 

Z B ( NET ) = EDS(NB) 

ZC(NE) = DV(NB) 

ZC(NET) = DVSAVE(NB) 

ZD(NE) = RHOEOP(NB) 

ZD ( NET ) = RHES ( NB ) 

ZE(NE) = RHOIOP(NB) 

ZE ( NET ) = RHIS(NB) 

600 CONTINUE 
PCI) = 5.0 
KKK(l) = 64 
KKK ( 2 ) = 2 
KKK ( 3 ) = NBB 
KKK ( 5 ) = NBB 
NB2 = 2*NBB 

CALL SCALE ( NB2 , Z A , KRSTR ) 

CALL PLOTMY (Z,ZA,KKK,P) 

WRITE (6,602) KRSTR 
ND = 0 

DO 606 NB = 1 , NBB 
NE = NBB - ND 
NET = 2*NBB - ND 
Z ( NE ) = XDV(NB) 

Z ( NET ) = XDS(NB) 

ND = ND + 1 
606 CONTINUE 

CALL SCALE ( NB2 , ZD, KRSTR ) 

CALL PLOTMY (Z,ZD,KKK,P) 

WRITE (6,603) KRSTR 
ND = 0 

DO 601 NB = 1 , NBB 
NE = NBB - ND 
NET = 2* NBB - ND 
ND = ND + 1 
Z ( NE ) = XDV(NB) 

Z(NET) = XDS(NB) 

601 CONTINUE 

CALL SCALE ( NB2 , Z E , KRSTR ) 

CALL PLOTMY (Z,ZE,KKK,P) 

WRITE (6,610) KRSTR 
ND = 0 

DO 611 NB = 1 , NBB 
NE = NBB - ND 
NET = 2* NBB - ND 
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ND = NO + 1 
Z(NE> = XDV(NB) 

Z ( NET ) = XOS(NB) 
fell CONTINUE 

CALL SCALE ( NB2 , Z B , KRSTR ) 

CALL PLOT MY (Z,ZB,KKK,P) 

WRITE (6,604) KRSTR 
ND = 0 

DO 608 NB = 1 , NBB 
NE = NBB - NO 
NET = 2*NBB - ND 
Z ( N E ) = XDV(NB) 

Z ( NET ) = XDS(NB) 

NO = NO + 1 
608 CONTINUE 

CALL SCALE ( NB2 , ZC , KR STR ) 

CALL PLOTMY (Z,ZC,KKK,P) 

WRITE (6,605) KRSTR 
IF ( IWRITE.EQ.l ) RETURN 
P(l) = NBB 
KODE = 64 
ND = 0 

DO 500 NB = 1 , NBB 
NE = NBB - NO 
Z ( NE ) = XDV(NB) 

ZA(NE) = RHODOP ( NB ) 

ZD(NE) = RHOEOP ( NB ) 

ZB(NE) = EOV(NB) 

ZC(NE) = DV(NB) 

ZE(NE) = RHOIOP(NB) 

ND = ND + 1 

500 CONTINUE 
WRITE (6,501) 

CALL SCALE ( NBB , ZA , KRSTR ) 

CALL PLOTXY ( Z , Z A , KODE , P ) 

WRITE (6,502) KRSTR 
WRITE (6,501) 

CALL SCALE ( NBB , ZD , KRSTR ) 

CALL PLOTXY(Z,ZD,KO()E,P) 

WRITE (6,505) KRSTR 
WRITE (6,501) 

CALL SCALE ( NBB , ZE , KRSTR ) 

CALL PLOTXY ( Z , Z E , KODE , P ) 

WRITE (6,510) KRSTR 
WRITE (6,501) 

CALL SCALE ( NBB , Z B , KRSTR ) 

CALL PLOTXY ( Z , Z B , KODE , P ) 

WRITE (6,503) KRSTR 
WRITE (6,501) 

CALL SCALE ( NBB , ZC , KRSTR ) 

CALL PLOTXY ( Z , ZC ,KODE , P ) 

WRITE (6,504) KRSTR 

501 FORMAT ( 2HPT ) 

502 FORMAT ( 2HPL , 47X , 8HNU ( X 10** , I 3 , 6H ) VS X) 

503 FORMAT ( 2HPL , 47 X , 7HE ( X 10**, I 3 , 6H ) VS X) 

504 FORMAT ( 2HPL ,47X , 8HDV ( X 10** , I 3 , 6H ) VS X) 

505 FORMAT ( 2HPL , 47 X , 8HNE ( X 1 0**, I 3 , 6H ) VS X) 

510 FORMAT ( 2HPL , 47X , 8HN I ( X 1 0** , I 3 , 6H ) VS X) 

602 FORMAT ( 2 HPL , 47X , 8HND ( X 1 0** , I 3 , 6H ) VS X/ 

12HPL,44X,20H+ R I CHARDSON-DUSHMAN/ 2HPL , 44X , 10H* SCHOTTKY) 
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603 FORMAT ( 2HPL * 47X * 8HNE ( X10** t 1 3 * 6H ) VS X/ 

12HPL 1 44X $ 20H+ R I CHARDS0N-DUSHMAN/2HPL »44X , 10H* SCHOTTKY) 

604 FORMAT ( 2HPL* 47X t 7HE ( X10** T I 3 t 6H) VS X/ 

12HPL»44X,20H+ R I CHARDS0N-DUSHMAN/2HPL ,44X , 10H* SCHOTTKY) 

605 FORMAT ( 2HPL 1 47X T 8HDV ( X10** f 1 3 , 6H ) VS X/ 

12HPL 1 44X * 20H+ R I CHARDSON-DUSHMAN/ 2HPL f 44X 1 10H* SCHOTTKY) 
610 FORMAT ( 2HPL t 47X, 8HNI < X10** f I3 T 6H) VS X/ 

12HPL,44X,20H+ RICHARDSON-DUSHMAN/2HPL,44X,IOH* SCHOTTKY) 
RETURN 
END 

$ I BFTC VINE 

SUBROUTINE VIN ( NV , DVV , T , ANS ) 

COMMON /MV/ P 1 1 EM» AK * AM 
C 

C IF NV - 1 COMPUTATION FOR ELECTRONS 

C IF NV = 2 COMPUTATIONS FOR ATOMS AND IONS 

C THE PROPER T AND DV IS SPECIFIED IN THE MAIN PROGRAM 

C 

IF (NV.EQ.l) GO TO 10 
C = SORT (2.0*AK*T/tPI*AM)) 

GO TO 20 

10 C = SORT (2.0*AK*T/(PI*EM)) 

20 XE =ASQRT(DVV/(AK*T) ) 

Y = ERF ( XE ) 

ERFC = 1.0- Y 

ANS = C*EXP(~ABS(DVV)/(AK*T) )/ERFC 

RETURN 

END 

$ I BFTC TRAPE 

SUBROUTINE TRAP ( X , DV ,NBB , ANS ) 

DIMENSION X(40)*DV(50),ANS<40) 

C 

C SUBROUTINE TO INTEGRATE TRAPEZOI DALLY 

C 

H = DV(2)/2.0 

SUME - 0.0 

DO 60 NB = 1 t NBB 

IF ( NB.EQ.l ) GO TO 50 

SUME = SUME+X ( NB— 1 ) +X ( NB ) 

ANS(NB) = SUME*H 
GO TO 60 

50 ANS(NB) = 0.0 
60 CONTINUE 
RETURN 
END 

SDATA 

1 

2 • 4E+3 
1 

3*0E+O 

1 

1 • E+l 3 
1 

2 • 5E+3 
1 

2.4E+3 
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APPENDIX C 


SYMBOLS FOR IBM OUTPUT SHEETS AND FORTRAN IV LISTING 

The symbols are presented here in the order of appearance on the output sheets. 


Output 

labels 

FORTRAN 

variables 

Symbols 

Description 

Units 

I 

AI 

I 

ionization poten- 
tial for plasma 
atoms 

V 

TE 

TE 

t e 

emitter temper- 
ature 

°K 

PHI 

PHI 

ecp 

work function 

V 

NEP 

EPN 

N 

ep 

plasma electron 
number density 

cm' 3 

TEP 

TEP 

T 

ep 

plasma electron 
temperature 

°K 

TIP 

TIP, TIPP 

T. 

iP 

plasma ion tem- 
perature 

°K 

LAMBDA 

AMDA 

X DT ep 

plasma Debye 
length 

cm 

PV 

PTEST 

p vp 

vapor pressure 
of plasma ele- 
ment at Tg 

torr 

(133. 322 
(N/m 3 )/torr) 

LAMBDA(TE) 

AMDATE 

x dt e 

emission Debye 
length 

cm 

DV 

DV 

AV 

sheath potential 
measured from 
plasma electron 
potential 

V 

ND(DV) 

RHODOP 

P AV 

net number density 
of charge at AV 

cm' 3 

NE(DV) 

RHOEOP 

P e 

electron number 
density at AV 

cm” 3 
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Output 

labels 

FORTRAN 

variables 

Symbols 

Description 

Units 

NI(DV) 

RHOIOP 

Pi 

ion number den- 
sity at AV 

cm’3 

E(DV) 

EDV 

e av 

electron electro- 
static field at AV 

V/cm 

X(DV) 

XDV 

X AV 

distance from emit- 
ter to AV 

cm 

JEE 

EJ 

ie 

emitted electron 
current density 

A/cm 2 

JEP 

EPJ 

^ep 

plasma electron 
random current 
density 

A/cm 2 

JIP 

PIJ 

•*ip 

plasma ion random 
current density 

A/ cm 2 

JAP 

APJ 

•*ap 

plasma atom equiv- 
alent random 
current density 

A/cm 2 

J 

CJ 

J 

net current density 
through sheath 

A/ cm 2 

PP 

PP 

P P 

plasma pressure 

torr 

(133.322 

(N/m 2 )/torr) 

JIE 

EIJ 

He 

emitted ion 

current density 

A/ cm 2 

JAE 

AEJ 

^aE 

emitted equivalent 
atom current 
density 

A/cm 2 

JA 

SJA 

^a 

net equivalent atom 
current density 

A/cm 2 

JI 

SJI 


net ion current 
density 

A/cm 2 

JE 

SJE 

U 

net electron 
current density 

A/cm 2 
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Output FORTRAN Symbols Description Units 

labels variables 


JA/JAP 

JE/JEP 

JI/JIP 

DVS 

XDVS 

NAP 

XD/LAM 

SC 

PHZ 


EDVS 

DVSRD 


DVS/RD 

ELM/RD 

PHZZ 


DRD/KT 

NTP 


SDJA 

SDJE 

SDJI 

DVS 

XDVS 

APN 

XLAM 

SC 

PHIZ 


EDVS 

DVSRD 


DVSPZ 

EDPZL 

PHIZZ 


DRDK 

TPN 


4V S 

X AV S 


(±0. 511x10® 
x 2. 82x10" 13 Eg) 1 / 2 

^o 


E AV g ~ E E 

<P -<P 0 


(c p nn = cp n for 
v ' oo ^ o 

equilibrium and 
electron sheath) 


•'a^ap 

^e^ep 

j^/j^p 

overall sheath V 

voltage drop 

effective sheath cm 

thickness 

plasma atom cm -3 

number density 

x s /x D 

Schottky depression V 

of work function 

plasma potential V 

(work function 
for no sheath) 

electrostatic field V/cm 

at emitter 

Richardson-Dushman V 

overall sheath 
voltage drop 

AVg/AVo = AVg /{(p - <p Q ) 

E E X D /(<p - cp Q ) 

plasma potential at V 

equilibrium (work 
function for no 
sheath and no net 
current) 

e\cp - cp Q \/KT e 

-3 

total particle num- cm 

ber density in 
plasma 
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Output 

labels 

NCE 

NTE 

RD/KTE 

X/LMTE 

ELT/RD 

NEPA 


NIPA 


FORTRAN Symbols Description Units 

variables 

-3 

CEN total charge number cm 

density at emitter 

-3 

TNE total particle num- cm 

ber density at 
emitter 

PHAT e\cp - <p 0 \/kT e 

XAMTE X s/ X DE 

EAMPHZ EgAj-jg/(cp - cp Q ) 

WNEPA NE(AV) at 


AV = 0. 0, ap- 
proximate value 
of N from 
sheath calcula- 
tions (Positive- 
Ion Sheath Pro- 
gram) 

WNIPA NI(AV) at 

AV = 0. 0, ap- 
proximate value 
of Nj from 
sheath calcula- 
tions (Electron 
Sheath Program) 

ND, NE, NI, E, 

DV, and X on 
plots correspond 
to ND(DV), 
NE(DV), NI(DV), 
E(DV), DV, and 
X(DV) in the pre- 
ceding list. 
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APPENDIX D 


IMPORTANT VARIABLES AND CONSTANTS IN FORTRAN IV LISTING 


The variables given in the following list are not included in the output. 


AK 

AM 

API 

C3 

C4 

CAT 

DVSAVE 

EDS 

EE 

EIJS 

EJS 

EM 

J1 

J2 

J3 

J4 

JA 

JB 

JC 

KOUNT 

PHD 

PI 

RHDS 


Boltzmann constant, k 
atom particle mass, m„ 

<x 

(Planck's constant)/(electronic charge), h/e 
gas constant 
electronic charge, e 

value of JA, JB, JC whichever is largest (Positive-Ion Sheath Program) 
Richardson-Dushman value of AV 
Richardson-Dushman value of E^y 
field at Schottky emitter, Eg or E^y 

S 

emitted ion current density j. g with Schottky correction (Electron Sheath 
Program) 

emitted electron current density j e g with Schottky correction (Positive-Ion 
Sheath Program) 

electron particle mass, m g 

emitted electron current density, j g g 

electron current from plasma that reaches emitter, j e p exp(-AVg//<T e p) 
ion current from emitter that reaches plasma, j. exp(-AVg//<Tg) 
plasma ion random current density, j. 

^kount " ^kount-1 

^kount ” ^kount-1 

^kount ” ^kount-1 
counter 

work function with or without Schottky correction, ecp or ecp ± e (±eEg) 

3. 14159 

Richardson-Dushman value of P^y 


1/2 
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RHES Richardson-Dushman value of p g 

RHIS Richardson-Dushman value of p^ 

RHOD (Net number density of charge) (Electronic charge) 

RHOE (Electron number density) (Electronic charge) 

RHOI (Ion number density) (Electronic charge) 

SE electronic charge 

TAP plasma atom temperature, T ap = T ip 

TENT 0. 1 percent of either Jl, J2, J3, J4, j , j whichever is smallest (Positive- 

a.p d© 

Ion Sheath Program); 0. 1 percent of either j g , j i ip , j ie , j , j ae which- 
ever is smallest (Electron Sheath Program) 
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Figure 2. - Flow diagram for 


























N* 

0 

1 


o 

O' 

o 


< 

a 

® 

x 

< 


U O 
I I 

UJ JJ 

in rg 
r- 
in o 
in 

® o 
»“* O' 


<r <r in in in ® in in in m in in rnn n o ^ 
oooooooooooo oooooo 


I I I I I I I I I I I I I I I I I I 

juuu-uojajajLU jLiuJajiXjvu uj jj uj lu iu UJ 

o O' O' inooHO'TOflM nmni njo® 

® oo ® ^ m n rf) n-o o' cr rg rn m ® in m ^ 


-r ® r» o or*- mo® ® ® rg *h m ® r-* o' o 
eg rg m <~h omi -oo h u in g'Mn O'UH 
o ® ® 'Tin'n>}-'rrMNf v H 0 'H- 0 'nno 
nrgo ^iNninmcoHino >r Oin h 


om>f niNHHNfn o 


rg ig ig ig nj in i\j i\) i\i ni m r) m m cn m m in m m 

O O U O O U O O O l_> Sj l_> u o o o u o o 


rg rg ® sj- 

000 o 

1 I I I 

aj UJ JJ UJ 

ooh m 

o O <-H ■—< 
ii'un ii\ 
<j<mn in 
OUOU'T NO 
l~ I- INJ — • 
• • • * 
mm r» ig 
I 


im p» ^ r» m ig t'i o' m 

r- m o' 'T m ui m ir in ni co r- >o in m 'T f w a >r 

i-» ^ p» 4 JUJI— i'-nu^ ujjj —• uj oj nt u> ^ r» u t 

O'HOP-Moniouir* h <y in r* n o p* ui 01 m 

m GO o r- m ® m jo -g — « — «oO U'O'O' OJ co aj ® 

m vr uw - 4 “ ni n> in in o h in m m >t m ® r- ~o O' 

I rH rN m vf *n ® r- oo o' H ^ H >— i ^ r-< h h h <-■« 


a®>» > 
«g. u 


o 

o 


m m m m cn rn rn m mfnenrnmrn<nrnfnmcnm m 


o O' 

o "• 

O ni 
O >fr 
O rr) 
O O 


O nj 
O' -a 
® o 

O -1 


« in g- 

4 O -I 
O nj 

g m 


rn 

O h- 
r- m 
m m 
rg o 

rg rg 


nT rn 
® r- 
® ® 
jo m 

o m 

m rn 


m rg 
o* o 
■O' O' 
O' r— 
O' 't' 
m vt 


® o 

4- O' 
^ 0 
N* JO 
O' 4- 

't* -n 


n m co -i h o -o 

co o n- o rn ■n jo 
^ g- h co g- h o' 
Hg>m Hnjin O' 

o m r- m o' m 
o-on s co® O' 


i i i t I 


I f i i I l i i l. I l I i I 


cgrgo 

??° 

JJ UJ uu 

in og 

OHN r- 

in O' vt 'H 
® m ® O' 
r- -n o rn 

sg co r- -» 
• • • • 

ni^rn'o 


o 

o 


eg 'g rg rg rg Mrgcgcg'grg-grgrgrgrgrgrgrgrg rg 


JJJJ JJ JJ UJ JJ UJ LU 

®®rgo'rgrgniO 

i ® in rg 


O' s|- ® JO O' 4- 

ins i — ~ ■ - 

«o <n ■ 

rg O' ® rn UJ l 


. O 

O' o rn 


O' o 

® m 

- 4 - O' 

*-l n}- 

rg ® 
® rn 


I UJ JJ UJ UJ l 


m rg i 
O' I s - r 
*H O v 

st- ® r 


> — ♦ rg 

i ® ® 

i 4 " m 
i in -i 
' in rg 


h co n n ~h ' 


'JJ LU 

® m 

Sl" rg 
m -h 

in rg c 
go <g r 

O ST r 


o 

rg in 


O' jo cyco-o r- ® ® ® nmn nmg g gg m 


— t rg O' r- 
ooo o 
I I I 

JJ JJ LU JJ 
MrgO h 
® - 4 " rg rg 
Oun rg 

g»«‘ n 

^ OJ IT( r -4 
rg J rg rn 


rgrgrgrgrgrgrgrgrgfMrgrgrgmrnmmmr 


1 JLlLUUJ UJLUUJ UJ JJUJ JJI 

>o n h *n to n o ornuir-i/'m<gr» 0 'r-o'r 


O' ® rg O' i 
® gr ig gj i 
grmiNnnHginingmo'iNi 


o h in in im >g g /I g g r x m in —• . 

I I I I I I I I I I I I I I 


i n hh rn 


i o ®r 
i m m s 


orgasm tMnHNfricgO'O'Nfflgg'inHN 
O' O' ® ® r* r* r» gj ® m u> ui g- g m n in >n in — t 
g "igr* g h mm in w< g m g n v h aj #i ig g« 
O*^i'g"grn'r>gin®'or»®o' 0 'o-H^igr»> n> 

oooo o ao ooooooo-h-h— »~<-h — < 
o o 3 0 0 OOQOOOOOO oooooo o 


JJ LU ill til 

h g rg o' 
o i— in gj 

ogrg H 

ig ® o o' 
m gj o' ig 


i »'i sr 
I I 


44 


NAP = i. 33 t SSCfc 13 XO/LAM = 1 .S 6542 JE UO P HZZ= 3.69222 EUVS = 1 . 9 E 3543 E C 3 

NT P = * . 33 1 599 E 13 NCE = 2 . 334936 E 13 ME = 4 . 671535 E 13 KL/nTE = 6 .i 2 S 372 E -01 

X/LMTE = i . 0 26364 E CC ELT/RO = 1 . 323154 E CO NtPA= 9 . 265926 E 12 


59€737£ C5 LAMBDA C T£ I = 1.0687E-04 


4 - 

o 

» 

UJ 

r- 

o 

O' 

o 



45 


(a) Numerical values. 

Figure 3. - Example output for Positive-Ion Sheath Program. 


1600. -1400. -1200. -1000. -800. -600. -4CC. -200. 0. 200. 400 



46 


ND(X10**-10i VS X 
81 CH ARUSuN-DUStiMAN 
SCHGTTKY 
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Nt ( X10** -9) VS X 
KI CH ARUS0N-UUSHM4N 
SCHOTT K V 

(b) Continued. Richardson-Dushman and Schottky results, 
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M (X10*«-10I VS X 
RI CH ARDS ON- DUS H MAN 
SCHUTTKY 
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■ RICHARCSuN-OUSHMAN 
> SCHGTTKY 

(b) Continued. Richardson-Dushman and Schottky results. 
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RICHARDSGN-DUSHHAN 

SCHOTTKY 


1600. -1400. -1200. -1000. -800. -6CC. -4CC. -2C0. 0. 200. 400 



(c) Schottky results. 
Figure 3. - Continued. 
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NE ( X10 ** -9 ) VS 


X*E 4 Y*E— 3 
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(c) Continued. Schottky results. 
Figure 3. - Continued. 
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(c) Concluded. Schottky results. 
Figure 3. - Concluded. 
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500. 1000. 1500. 2000. 2500. 3000. 3500. 4000. 4500. 5000 
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M>mo**-m vs x 

RICHAROSON-QUSHMAN 

SCHOTTKY 
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NE ( X10 **—11 ) VS X , * . . _ . , _ .... .. 

h i cHAHosuN-uus hm^ n (b) Richa rdson- Dush ma n and Schottky results. 

SCHOTTKV 

Figure 4 . - Continued. 
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NASA-Langley, 1968 
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DV1X10** o) VS X (c) Concluded. Schottky results. 

Figure 4. - Concluded. 



